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Abstract 

The conceptual difference between equilibrium and non-equilibrium steady 
state (NESS) is well established in physics and chemistry. This distinction, how- 
ever, is not widely appreciated in dynamical descriptions of biological populations 
in terms of differential equations in which fixed point, steady state, and equi- 
librium are all synonymous. We study NESS in a stochastic SIS (susceptible- 
infectious-susceptible) system with heterogeneous individuals in their contact be- 
havior represented in terms of subgroups. In the infinite population limit, the 
stochastic dynamics yields a system of deterministic evolution equations for pop- 
ulation densities; and for very large but finite system a diffusion process is ob- 
tained. We report the emergence of a circular dynamics in the diffusion process, 
with an intrinsic frequency, near the endemic steady state. The endemic steady 
state is represented by a stable node in the deterministic dynamics; As a NESS 
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phenomenon, the circular motion is caused by the intrinsic heterogeneity within 
the subgroups, leading to a broken symmetry and time irreversibility. 

Keywords: circular motion; multi-dimensional birth-death process; non-equilibrium 
steady state; Omstein-Uhlenbeck process; time irreversibility. 

1 Introduction 

The mathematical description of population dynamics is universally based on 
"population change per unit time = birth rate — death rate." [[Tl |2]|. One of 
the most widely employed nonlinear models for teaching population dynamics 
is n = rn{l — n/C), where r is the growth rate per capita when the population 
is infinitely spars, and C represents carrying capacity. However, a careful reflec- 
tion immediately leads to the following question: Does r(l — n/C) represent a 
decreasing per capita birth rate; or is the birth rate being linear rn and death rate 
being rn^ / Kl Such an explicit distinction has no consequence in ordinary differ- 
ential equation models. However, it has important consequences to stochastically 
fluctuating dynamics of finite populations [3 J. 

The distinction is necessary when one develops an individual behavior based 
stochastic model. Population dynamics in terms of stochastic, distributed individ- 
ual behavior is a more refined description of the biological reality. The nonlinear 
dynamics one often observes is a collective, emergent behavior at the level of a 
population within which individuals making seemingly random choices with bias. 
In sociology and economics, this type of modeling is called agent-based, and in 
finance it is called behavior finance. 

A case in point is the infectious epidemics among a large population with 
several behavioral subgroups. It has long been recognized that various types of 
heterogeneities play fundamentally different roles in epidemics such as sexually 
transmitted diseases (STDs) [4]. Statistical and analytical studies have defined 
and characterized heterogenous phenomena such as the existence of core groups 
leading to the "80/20 rule" [Sll. 



2 



The rise of "network theory" and its graph representations have provided a 
new perspective on modeling heterogeneity in population dynamics. In particular, 
networks with scale-free characteristics have been used to represent diverse popu- 
lation interactions, such as human sexual contacts air transportations fTTll, 
and the World Wide Web |l9l[T0l|. Dynamics of such networks has attracted much 
attention flUIIlIIlllIl. 

Another widely used approach to heterogeneity is to introduce a probability 
distribution over a population. This naturally leads to dynamic models based on 
stochastic processes. The stochastic approach has received growing attentions in 
recent years. For example, Nasell has extended the classical methods of moment- 
closure and pairwise approximation (see lfT6l and the references cited therein.) 
Stochastic approaches have also been combined with network dynamics [fTTl [TTl 

In the present work, we consider the statistical behavior, at the individuals 
level, which can itself be "heterogeneous". In other words, each subgroup has a 
different "contact number". This notion is motivated, though no need to be, by the 
degree distribution in the network theory. Specifically, using nodes to represent 
individuals and edges to represent contacts between individuals, we envision every 
individual has a fixed number of half-edges representing his/her level of social 
activity. The heterogeneity in our model, thus, comes from the different number 
of half-edges an individual possesses. A contact between two individuals is then 
represented by a connection of two half-edges. 

Compared with relatively long time period of being infectious (/) and suscep- 
tible {S), the act of disease transmission is often very short, and will be considered 
as instantaneous. Therefore, we shall assume that the contact between two half- 
edges as a collision event between two individuals. Such contacts can lead to an 
individual changing state from 5 to /. In this way, our model is different from 
the pairwise models which deal with concurrent partnerships [fTTl [TTl [T8l[T9ll20l l. 
There exists at most one "collision" in the system at any time. 

These considerations naturally lead to the mass-action law, which we adopt: 
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With an identical rate, any two half-edges from different individuals are possible 
to collide with equal probability. The actual collision events are stochastic, hap- 
pening in serial and independently. This is known as a Poisson flow. Furthermore, 
we assume that the recovering time of an infectious individual obeys an exponen- 
tial distribution. Based on these assumptions, we define a stochastic process on 
the level of population subgroups. For a system with K subgroups, the stochastic 
process is a birth-death process with multi-dimensional discrete state space 

m. 

In the limit of infinite population size, using the method of ^7-expansion (22^, 
we show that our stochastic model becomes a system of deterministic evolu- 
tion equations which are closely related to the models previously studied by Laj- 
manovich and Yorke flUl, May and Lloyd lfT3l . and others [|T4ll . This indicates that 
the assumptions we made for the mechanism at the individual level are reason- 
able; and in fact we have obtained a stochastic counterpart of the classical model. 
A stochastic model requires more explicit assumptions than those for differential 
equations. 

As for the cause(s) of often observed noisy oscillations (fluctuations) in epi- 
demiological data, it is still controversial . Deterministic framework focuses 
on interactions between external forcing and inherent frequencies in nonlinear dy- 
namics [12411 . Stochastic models, however, illustrate a fundamental role of intrinsic 
randomness in the patterns of disease cycles [|25l |26l |27ll . Such a debate echoes 
the nature of fluctuations/oscillations in the concentrations of chemical species in 
a small volume. It has been shown that in a stochastic, diffusion process model 
one can unambiguously distinguish two types of mechanisms [28 1: stochastic fluc- 
tuations and nonlinear dynamic complexity. Furthermore, it has been shown that 
nonlinear oscillations in a stationary stochastic dynamics are intimately related 
to a probability circulations associated with time-irreversible Markov processes 
[|29l[30ll3T| . One of the objectives of the present work, thus, is to call for atten- 
tions to population dynamic studies of chemical species which might serve as a 
paradigm for dynamics of more complex systems [[3l[32|[. 



4 



The interplay between stochasticity and nonlinearity in epidemics was stud- 
ied in [|26ll : an oscillatory spiral type steady state in the deterministic system was 
shown to be amplified by the demographic stochasticity. In our system, the oscil- 
latory dynamics is a consequence of finite population; it disappears in the deter- 
ministic nonlinear dynamics. 

Classical SIS systems with well-mixed homogenous individuals have no oscil- 
latory dynamics in either deterministic or stochastic models. In the former, there 
exist at most two attractors: a trivial stable node and a non-trivial stable node. The 
term "node" here refers to a type of fixed point with real eigenvalues in dynamical 
systems theory; it should not to be confused with the same term in graph theory. 
(A planar fixed point with complex eigenvalues is called a "spiral".) Stochastic 
SIS model is a one-dimensional birth-death process [33J; it is known that such a 
process is time-reversible if it is stationary. 

Deterministic SIS dynamics with heterogenous individuals is multi-dimensional; 
still it exhibits the same type of behavior as the homogeneous case flU [141: All 
fixed points are nodes. However, for the corresponding stochastic model in this pa- 
per, we discover that the multi-dimensional birth-death process can exhibit circular 
dynamics. We shall first demonstrate this by a linear theory near the non-trivial 
stable node, i.e., via an Ornstein-Uhlenbeck (OU) process. We then investigate the 
nonlinear, multi-dimensional birth-death process and show the circular motion as 
a fundamental property of the stochastic epidemic process. It is a consequence of 
the heterogeneity among individuals. 

The paper is organized as follows. Sec. 2 describes the stochastic contact 
process among individuals leading to the definition of a multi-dimensional birth- 
death process. In Sec. 3, we carry out the ^7 expansion of van Kampen [|22ll and 
obtain a diffusion approximation of the multi-dimensional birth-death process in 
the limit of large population size N. The Q expansion consists of two steps which 
are generalizations of the Law of Large Numbers and the Central Limit Theorem 
[|34l . In the first step a system of ordinary differential equations (ODE) is obtained 
under the A^^^ scaling. Then conditioned on an ODE solution and with a N~^^'^ 
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scaling, an OU process is obtained [|3l[37l. 

In Sec. 4, the deterministic ODE system is analyzed. We show under different 
conditions the system yields either a positive non-trivial stable node or a trivial 
one at the origin. The former corresponds to an endemic state while the latter 
corresponds the infection being eradicated. 

Then in Sec. 5, the properties of the OU process near the positive non-trivial 
stable node are studied. While the stable node gives no indication of any oscil- 
lation by the ODE, the OU process undergoes a sustained circular motion around 
the stable node — It is represented by a nonzero stationary circular flux. 

In Sec. 6, we use the simple system of two subgroups to further illustrate the 
circular pattern of infectious dynamics. We show it is not a result of fi-expansion 
approximation; we demonstrate that the discrete multi-dimensional birth-death 
process violates the so-called Kolmogorov cycle criterion. Therefore, according to 
a mathematical theorem in irreversible Markov processes, the circular flux exists 
in the discrete model — This feature originates from the heterogeneous contacts 
and it does not occur in SIS systems with homogenous populations. The paper 
concludes with Sec. 7. 

2 The model 

2.1 Individual contacts and recovery 

We consider total individuals in a population. We assume that each individual 
has a fixed "contact degree" in an epidemic, which is represented by the number 
of half-edges (or valency) associated to the individual. We further assume that 
the entire population can be partitioned into subgroups according to the contact 
degrees of individuals. For example, groups 1, 2, 3, etc., represent individuals 
with contact degrees between — 10, 11 — 20, 21 — 30, etc., respectively. The 
rate of encounter between two individuals from subgroups i and j is assumed to 
be proportional to i x j. 

Let Nk be the number of individuals in the /c*^ subgroup and total population 
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= "^j^^i N^. We shall denote the fraction of the population N^/N by D^; 
k = 1,2, ■ ■ ■ , K, with K being finite. In graph theoretical language, the Dk is the 
degree distribution for the contact network: XIa^i -^fc = 1- 

Following the standard compartmental modeling of epidemics, we assume ev- 
ery individual in the population is in one of the two states: susceptible (denoted by 
S) or infected (I). In the present work, an individual becomes infectious immedi- 
ately after being infected. Let Tj denotes the infectious period of individual i, and 
we assume that all T^; i = 1, 2, ■ ■ ■ , A^, are independent and identically distributed 
following an exponential dwell time with mean value 1. An infectious individual 
becomes 'susceptible' once again when his/her infectious period is terminated. 
These assumptions imply the heterogeneity we consider is in the contacts between 
individuals, while the infected individuals behave statistically homogeneous. 

One can visualize the contacts among individuals in our model as follows: 
There are regular, repeated "touching" or "collisions" within pairs of half-edges, 
which trigger one of the individuals with certain probability to change its state. Ev- 
ery "touching" is instantaneous. If one of the individual is infectious and the other 
is susceptible, then an infection can occur. In some sense, the dynamics is no dif- 
ferent from an autocatalytic chemical reaction in aqueous solution A + X 2X, 
as in the celebrated Lotka-Volterra model ll35l . See |l3l for a discussion on dy- 
namic isomorphisms between cellular biochemical and epidemiological popula- 
tion dynamics. 

2.2 Defining a stocliastic epidemics 

We now turn the previous verbal description into a mathematical model. We de- 
note the stochastic demographic process 

X(t)= - - t >0, (2.1) 

in which I^^^ (t) is the number of infected individuals in the k^'^ subgroup at time t, 
and accordingly Nk — /^'^•'(t) will be the number of susceptible individuals. X{t) 
is a continuous-time, i^'-dimensional birth-death process. Since < I^''\t) < N^, 
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X{t) takes values in a bounded subset of the X-dimensional lattice Z^. 

When initially there are infectious individuals in a population, there will be 
a stochastic epidemic process. Let A be the rate of infection for a single half- 
edge. A/(^^^^ rnNm) then is normalized by the total number of half-edges in 
the system such that in the infinite population limit, Nm, A ^ oo but the ratio is 
unchanged. Then, a give half-edge makes pair with an infectious one at the rate of 

For a given susceptible individual with degree k, it will be 'infected' at the rate 
of k\. Thus the number of susceptible individuals — I'^^^t) decreases by 1 

in the fc*'' subgroup with rate k{Nk — I^''^)X. Illustratively, the transitions and 
corresponding rates of the process X{t) are given as follows: 

transitions: ± rates 



(2.3) 

K 

-ek = ( p,-- - ,0,-1 , 0, • • • , 0) = 

k 

Here k = 1,2,-- - ,K. Recall that without loss of generality the recovery rate of 
an infected individual is 1. 

Denotes the probability P{p;t) = Pr{X{t) ^ p}, p e Z^. The stochas- 
tic dynamical system is characterized by the master equation for the probability 
distribution 

d ^ 

-P(0,t) = J2j1{6+ek)P{0 + ek), (2.4a) 



fe=i 

K K 



jP{p.t) = J]4(p-e,)P(p-efc) + 5^J!(p + e,)P(p + e,) 



fe=i fc=i 

K 



-J2{J-iP) + JliP))PiP)^ {0<P<Nk) (2.4b) 



k=l 

K K 



'^P{Nk,t) = Y.JliNk-ek)P{Nk-ek)-J2j-(^k)P{Nk). 

k=l k=l (2.4c) 



dt 
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It should be noted that the solution to Eq. (12.41 ) with initial data P(p; 0) = 5{p — 
pi) is the transition probability Pr{p(t) = p|p(0) = pi} for the Markov process 
X{t). 

It is obvious that the is an absorbing state. This means the infinite long-time 
behavior of the system is always the extinction of the infectious population [i36ll . 
We are, however, interested in the pre-extinction dynamics. In many biological 
population dynamics, extinction is too long a time scale to be relevant. To explore 
the pre-extinction dynamics, we shall first consider large population {Nk) limit 
and establish the corresponding nonlinear deterministic behavior in the following 
two sections. Pre-extinction dynamics can also be studied in terms of a transient 
landscape [ITTl . However, we shall not pursue this approach in the present paper. 



3 Population-size (Q) expansion for large 

With increase population size in an increasing space, the population dynamics 
for discrete individuals can be described by a continuous variable representing 
the "density", or "frequency", or any other intensive quantity. The Law of Large 
Numbers from the theory probability then dictates the dynamics of the system ap- 
proaches to a deterministic behavior. Mathematically, a systematic method called 
VL expansion is widely used in statistical physics U2\ . See its recent application in 
epidemiological modeling [|26l . 

Let us first introduce the notation of a "step operator" [|22l . which repre- 
sents the elementary changes of the discrete X{t). It is defined by its transforma- 
tion of a function /(p): 

Efc/(p) = /(p + e,) and V(p) ^ f{p- e,). (3.1) 

Then the master equation in Eq. (12.41 ) can be rewritten in a much more compact 
form 

, K K 

-p{p- t) = j2 (E. - 1) j'-imp) + E (^^^' - 1) 4(p)^(p1- (3-2) 
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Now we introduce a new vector y, which scales as ^ so represents the density 
vector of infected population, and a fluctuation ^ — which scales as N^^/'^ such 
that 

Pk = Nvk + N'/^iu + o . (3.3) 

Accordingly the distribution is now written as a function of ^, 

P{p-t)=Ii{lt), (3.4) 
and the master equation (12.41 ) in the new variable takes the form 



dt ^ dt dSk 

k 

K 



Kvk + N-y^ii) 



ik) 

K 



m) 



+ E (^'^^'1^ + + ■ ■ ■) i^y^ + ^^''^^) n(a, (3.5) 

in which parameter 

w=^l^=i:*^*- (3.6) 

The conditions for the terms of order A^^/^ to vanish are 

^^ = -yk{t) + ^ADk-yk{t))Y,3y3it), k = l,2,---,K. (3.7) 

j=i 

This is a system of A'-coupled nonlinear ordinary differential equations (ODEs) 
for ykit): The deterministic nonlinear dynamics for infinite size population. The 
system of ODEs is a generalization of the one-dimensional logistic equation y = 
—y + \{D — y)y, which has two fixed points, one at and another positive one 
at D — 1/A when \D > 1. This result should be understood as the Law of Large 
Numbers for the stochastic process X{t). 
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The terms of order in (13.51 ) yield 

dt ^\ {k) 

This is a time-inhomogeneous linearly diffusion process centered at the time- 
dependent y. It should be understood as the counterpart of the central limit theo- 
rem for X{t). 

Eqs. (13.31 ). (13.71 ) and (13.81 ) together show that the evolution of the stochastic 
dynamical system described by master equation (13.21) . when the population size 
is large but not infinite, can be characterized by two parts: One represents a de- 
terministic nonlinear dynamics of the densities, yk{t), in infinite population-size 
limit, and another represents the fluctuations in terms of stochastic diffusion, with 
density function n(^fc, t), centered around the deterministic trajectory. It is impor- 
tant to point out that the yk{t) is the behavior of an infinitely larger population; 
it is not in general the mean dynamics for a finite population with nonlinear in- 
teraction. Mean dynamics of a finite population is of course "deterministic", but 
usually it can not be described by a self-contained set of autonomous nonlinear 
differential equations. See [[38l for a recent study on the moment closure problem. 




4 Deterministic dynamics with a stable node 

We now study the system of ODEs (13.71) . It should be noted that similar equations 
have been proposed phenomenologically in the past flU [121 [Ml. In particular, 
the elegant mathematical analysis in [j4l has dealt with a more general class of 
problems, which can easily be applied here. In order to make the material more 
accessible to readers with less mathematical backgrounds, here we present a sim- 
plified recapitulation of the earlier work for the particular ODE system (13.71 ). Our 
method is also more explicit; it can be conveniently adopted in numerical com- 
putations for the non-trivial stable steady state: Using the C and G{C) defined in 
Eqs. (14.11 ) and (14.21 ). one can easily obtain C* from arbitrary chosen initial C, with 

11 



computational iterations with rapid convergence. For more rigorous mathematical 
treatment, however, one is referred to flU]. 
Let 

C(t) ^ ^"^=^,!^^"^'^\ (4.1) 
{k} 

where < C(t) < 1 since mDm = (k) and Dm > Vmit). Um is the density 

of infectious individuals with degree m, so C{t) represents the mean fraction of 
infectious half-edges in the system. 

According to Eq. (13.71) . the positive equilibrium is located at yl = x+Afcc* ' ^ ~ 
1,2,--- ,K. Ifykit) < fgig, ykit) will increase, so does 

Even though the dynamics is if-dimensional, it is easier to work with the 
single variable C(t). We introduce a function of C: 

G\C) = TTT / ; — t;- (4-2) 

It can be explained as the averaged infectious half-edges caused by C in an unit 
time. We see that C (t) increases if 

C{t) = ^^=i7"^"-(^) < G{C). (4.3) 

{k) 

Alternatively, C{t) decreases when C > G{C). The auxiliary function G{C), 
therefore, has a fixed point C* 7^ as the solution of G{C) = C. Accordingly the 
steady states of Eq. (13771) . yl = y^fffr, can be obtained. 

The fixed point C* is determined by the properties of the function G{C). We 
note that G'(O) = 0, and 

Furthermore, we have G{C)'s derivative: 

(k) (1 + XmCy 

and G"{C) < 0. All together, we see that function G{C) is monotonically in- 
creasing and concave in the region < C < 1. Fig. [T] shows the function G{C). 
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There are two cases according to the values of A. When \ < \c = {k) / (fc^), 
G'{0) < 1. Then there is only one fixed point C* = and correspondingly yl = 
\/k. It is stable. In this case, the disease will fade out eventually. 

If A > Ac, G"(0) > 1, the fixed point becomes unstable. Simultaneously there 
emerges a unique stable fixed point < C* < 1, and accordingly yl = Y^§^§r. 
In this case, there will be an endemic state. The results can be seen directly from 
Fig. [B There is a trans-critical bifurcation at A = Ac, as in the logistic equation 
y = —y + X{D — y)y when A = 1/D. 




Figure 1: The functions G{C), given in Eq. (14.21 ). with two different values of A. When 
A > Ac, there will be a positive, stable node corresponding to an endemic steady state 
with a sustained infection. 



4.1 Linear analysis of steady states 

One can linearize equation (13.71) near its steady states and obtain 

f^BS (4.6) 
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in which the linear community matrix B has elements 



(k) - \- {k) 



+ 1 . 



(4.7) 



At the non-trivial, internal steady state, one can further obtain the complete 
eigenvalues of B: 



Ai = -AC*, Afc = -\kC* - 1; 2<k<K. 



(4.8) 



Therefore this steady state is a locally stable node. This means that there is no any 
indication of oscillatory behavior, at this deterministic level, near the steady state. 

At the steady state for extinction: y* = 0, we have B^i = XklDk/{k) — 5ki 
from (14.71) , which yields the eigenvalues 

.2 



, Aj = -1; 2<i<K. 



The eigenvector corresponding to the first eigenvalue Ai is 



K 



' D,{k^) 



+ A > 0, 



> . 



i=2.3,---,K 



(4.9) 



(4.10) 



In the neighborhood of 0, this eigenvector locates within the first quadrant. It is 
stable when A(A;^) / (k) <l, and unstable when A(A;^) / (k) > 1. 



5 Non-equilibrium fluctuations near the positive sta- 
ble node 

We now turn our attention to stochastic dynamics near the stable internal steady 
state y = y*. Note that we have shown through the above linear analysis that y* 
is a node, with all eigenvalues being real. To consider stochastic fluctuations, we 
consider Eq. (|3.8I) . Substituting the deterministic steady state solution into Eq. 
(13.81 ), it is reduced to a time-homogeneous Fokker-Planck equation defined on the 
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entire 



du{^, t) ^^n(e, t) - i?f n(e, t) ] , (5.1) 



dt \2 

whose stationary solution is called Ornstein-Uhlenbeck (OU) process. 

For the trivial stable node, the corresponding Fokker-Planck equation is only 
defined on the first quadrant. This leads to a singular diffusion equation, which is 
outside the scope of the present paper and will be the subject of a separated study. 

The constant drift matrix in Eq. (|5.1I) is precisely the community matrix B in 
Eq. (14.61 ). The diffusion tensor is a diagonal matrix Ak^iv*) = 2?/^ since we have 
the useful relationship ^{Dk — Vk) ^JUj = Ul from Eq. (13.71 ). The stochastic 
trajectories of OU process defined by Eq. (|5.1I) follows a linear stochastic differ- 
ential equation: 

^ = i?e(t) + rc(t), (5.2) 

where TV'^ = A, C(t) is a Gaussian white noise with E[C{t)C'^{t')] = 5{t - t')I, 
and / is the identity matrix. 

If one sets the initial condition 

K 

n(e,o) = n^te-cn, (5.3) 

i=l 

Then the fundamental solution, i.e.. Green's function, to Eq. (15.11) is a multivariate 
Gaussian distribution ||39ll40l. 



n(e;t) =noexp|-i (f-/2(t))^S(t)-i (j-f2(t)^Y (5.4) 
where Ho is a normalization factor. The first and second order moments satisfy 

'-f^-j:Bk,,k, '^^BS^EB^ + A. (5.5) 

j 

The solutions to Eq. (|5.5I) are: 

m = e*^eo, m = re(*-*')^Ae(*-*')^"rft'. (5.6) 

Jo 
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Since the eigenvalues of B are all negative, we have (^)(oo) = and ap- 
proaches to a as the solution to Eq. (|5.7I) : 

BE' + E'B'^ = -A. (5.7) 

Equation (15.71) is analogue to the well-known Einstein's fluctuation-dissipation 
relation, in which A is the covariance of the fluctuating white noise, B is the 
dissipative linear relaxation rates, and E' is the equilibrium covariance. In thermal 
physics, is proportional to ksT. 

Very interestingly, we note that A^^B is not symmetric. This implies the 
stationary process has certain breaking symmetry with respect to time reversal 
[|40l . In physics, there is a refined distinction between an equilibrium and a non- 
equilibrium steady state. Such a distinction is not widely appreciated in dynamical 
descriptions of biological populations in terms of differential equations in which 
fixed point, steady state, and equilibrium are all synonymous. 

According to the method previously developed in physics [|3l |40l SH, there 
exists a stationary circular current density 

= ($0^(0, (5.8) 

in which is a K x K matrix with zero trace and j^^^ is a divergence free vec- 
tor field V ■ j^'^-' = 0. The linear force B^ near the steady state, then, can be 
decomposed into two orthogonal parts 

B = --AE-^ + (5.9) 
2 

In Eq. (15.91) . the first term — ^AE~^^ is a conservative force with a potential func- 
tion f/(0 = I'C^ i'^'^y^ ^- It generates a motion towards the origin. In another 
words, it represents a stability landscape U{^). The second term generates 
a stationary circular motion on the level set of U. To see this, we observe the 
orthogonality 

jW±vn(0 (5.10) 



16 



due to = 0. The stationary probability distribution does not change 

along the direction of The maintenance of the stationary distribution comes 
from two parts: one is the detailed along the gradient of U{^), and another is 
circular along its level curve. 

For heterogeneous SIS with two subgroups, i.e., K = 2 and k = 1, 2, the 
matrix B from Eq. (|4.7I) is 



B 



and matrix 



/ X{Di-2yl-2y*2) . 2X{D-,-y*) \ 

D1+2D2 D1+2D2 

2\{D2-yl) 2\{2D2-yl-Ay*2) 1 

\ D1+2D2 D1+2D2 / 



(5.11) 



2y\ 

^=1 I, (5.12) 

2yl 

in which 7/2 ^ (0, D2] is a root of 

4 {ylf - (^2D, + I2D2 - y; + 20^ {d, + AD^ - = 0, (5.13) 

and 

yl = -^§^^-2yl. (5.14) 
2\[D2 - y2) 

Therefore, with parameters A = 1 and Di = D2 = \, thus {k) = 1.5, the matrix 

/ -0.0644 0.1204 \ 
$= . (5.15) 

y -0.1903 0.0644 J 

The eigenvalues of this $ is pair of pure imaginary numbers ±0.137i. This reflects 
an intrinsic frequency of the dynamics. Fig. |2K shows how the frequency changes 
with (k) for several different values of A. For heterogeneous populations with two 
subgroups A^^i and N2, the maximal effect of heterogeneity has to be (k) between 
1 and 2; where the imaginary eigenvalues of $ is also the largest. 

One of the widely used technique to capture and characterize oscillatory dy- 
namics is the method of power spectrum [[30|. For a strong oscillatory motion with 
noise, its power spectrum exhibits an off-zero peak at the intrinsic frequency. By 
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Figure 2: (A) For K = 2, the intrinsic frequence, i.e., the pure imaginary eigenvalues of 
as a function of ik) and A. {k) changes from 1.999 to 1.001 due to Di changes from 
1/1000 to 999/1000 , and the corresponding D2 = I - Di changes from 999/1000 to 
1/1000. (B): the IR (imaginary-to-real) ratio, the ratio of the imaginary eigenvalues of 
<^ to the sums of two real eigenvalues of matrix B, as a function of {k) and A. Note the 
critical Ac = 1 here; hence when A < 1 the endemic steady state disappears. 

"strong", we mean the intrinsic frequency has to be significantly greater than the 
relaxation rate of the stochastic dynamics. Il42l has shown that the ratio of imagi- 
nary part to real part of an eigenvalue has to be at least greater than = 0.577. 
Fig. shows the ratio between the intrinsic frequency, given by the imaginary 
eigenvalue of $, and the relaxation rate, given by the eigenvalues of the commu- 
nity matrix B. This explains why one does not observe an off- zero peak in the 
power spectra of the simulated dynamics (data not shown). 

6 Intrinsic circular dynamics in multi-dimensional 
birth-death processes 

We now address a crucial question unanswered in Sec. 5: Whether the novel 
circular dynamics near the endemic steady state, shown by the OU approximation, 
is a consequence of the population- size expansion approximation, or is it a general 
result for the finite population birth-death process. The answer is affirmative. To 
illustrate this, we shall again consider the special case oi K = 2 and analyze 
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the planar system in some detail. The same methods and results can easily be 
generalized to any K > 2, but the algebra will be more cumbersome. 

First, we need to address the issue of long-time behavior of the original birth- 
death process. In Sec. 4 we have shown that when A > Ac, the deterministic 
dynamics has a stable internal node representing an endemic steady state, while 
the extinction state Hi = Q {I < i < K) h unstable. On the other hand, it had been 
shown in Sec. 2 that the infinite long-time behavior of the original birth-death 
process is always extinction, even though it might take extremely long time. This 
disparity between the stable deterministic steady state and the stochastic long- 
time behavior indicates a separation of two different time scales: The time scale 
on which the infectious dynamics reaches a stationary pattern among different 
populations, and the time scale for extinction of the infectious population [f36ll43l 
l44l . In the light of the time-scale separation, this is a well-understood subject 
in population dynamics. There are in fact two fundamentally different types of 
extinction: one that occurs as a consequence of nonlinear dynamic, and another 
that occurs as rare events that is impossible according to nonlinear dynamics [[37l . 

To study the "long-time" dynamics in the pre-extinction phase, one of the 
methods of attack is quasi-stationary approximation [|36l [371 |43l . Conceptually, 
we only consider the stationary probability distribution conditioned on the sur- 
vival probability at time t. To carry out the computation, one can introduce a very 
tiny transition probability e ^ 1 representing the process can "return from" the 
state 0. This approximation abolishes the absorbing state. Since the total num- 
ber of states is finite for our original birth-death process and it is irreducible, an 
unique stationary distribution exits. Biologically, the e can be interpreted as an 
infinitesimal invasion by migration, or recurrence, of infectious individuals. The 
time evolution of probability in (|2.4h ) then becomes 

d ^ 

-P(0, t) = Y, ^-(0 + efe)P(0 + e,) - £5(P(0, t) e [I - e, 1]). (6.1) 

fc=i 

In the case of = 2, the probability transition rates in the planar system 
are illustrated in Fig. [3l Such a diagram is known as a master equation graph in 
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Figure 3: The "phase plane" for a two-dimensional birth-death process (/i,/2) where 
Ii and I2 take non-negative integers. All the possible transitions among four neighboring 
states are shown as the arrows and the corresponding rates are labeled by the transitions 
{i.e.,Eq. Q). 

stochastic chemical kinetics [|45ll46l . It is useful tool in "visualizing" the dynamics 
of the master equation in Eq. (I2.4h . on a par with the phase portrait of planar 
nonlinear dynamics. 

According to a key theorem in the theory of irreversible Markov chain, a suf- 
ficient and necessary condition for the existence of circular flux in the stationary 
process is the Kolmogorov cycle condition ||29l[30ll . In Fig. [3]the values for all the 
transition rates around a square cycle are given. Let denote the "clockwise" 

circular rate product 

j{m,n) 

<-'+ ~ Q{m,n)—^{m,n+l)Q{m,n+l)'^{m+l,n+l)Q{m+l,n+l)-^{m+l,n)Q{rn+l,n)~>-{rn,n) j 

(6.2) 

in which the q(m,n)^(m,n+i) is the transition rate from grid point (m, n) to grid 
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point (m, n + 1). Similarly J^™'"^ is the "counterclockwise" rate product 

j{m,n) 

<-'— ~ 9(m,n+l)— >-(m,n)'?(m+l,n+l)— s>(m,n+l)'?(m+l,n)^(m+l,n+l)'?(m,n)— s>(rn+l,n) • 

(6.3) 

Then the Kolmogorov cycle criterion 6 for each little square cycle is: 

j{rn,n) 

Substituting the values for the transition rates in Fig. [3l we have 6'(m „) = ^^^n+i ' 
shown in Fig. IH We note that 6{rn,n) is very close to 1 for large m or n. So 
generally the circular motion is weak. This explain why there is no circular motion 
in the deterministic dynamics of densities, when N ^ oo. 

The 9(jri,n)^ in Fig. |4]can be heuristically understood as "vortex" in a vector 
field. They are the fundamental units of circular dynamics. For a larger closed cir- 
cular path, its is simply the product of all the ^^'s within. The significance of the 6 
value lies in the theory of stationary, irreversible Markov processes. According to 
the cycle decomposition theorem for irreversible Markov processes [|29ll , the 0{rn,n) 
is the ratio of numbers of occurrence, following the stationary, stochastic birth- 
death process, of the clockwise cycle (m, n) — > (m, + 1) — )■ (m + 1, n + 1) — )■ 
(m + —7- (m,n) to that of the corresponding counter-clockwise cycle. A 
stationary process is time-reversible if and only if each and every cycle has 9 = 1. 

The Fig. |4]gives the insight that the infection dynamics is "moving clockwise" 
in the phase plane in its stationary state. We know that a one-dimensional station- 
ary birth-death process is always time-reversible. Therefore, irreversibility in the 
present example comes from the heterogeneity in the infectious population. This is 
a novel dynamical behavior emerged in population with heterogeneous structures. 

In order to see things more clearly, we carried out stochastic simulations for 
the stochastic process. Firstly we consider two subgroups with k = 1, 2. Other 
parameters used are A^^i = 50, A^2 = 50 and the transmission rate A = 1. The 
initial conditions are chosen randomly. Fig. [5] shows the typical trajectories of the 
two infecting population sizes, Ji and I2. For comparison, the two thick lines are 
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Figure 4: Stationary birth-deatli process can have a hidden circular dynamics, analogous 
to the vortex in a fluids |i47']. The existence of stationary circulation can be determined by 
the Kolmogorov cycle number ^(m,n) ■ For our model with K = 2, the 0(m,n) = m+2ra+i 
and the circulation is always clockwise. For a larger, arbitrary closed circular path, the 
value of is the product of all the individual 6''s within. A stationary Markov process is 
time-reversible if and only if all the 0's are 1. In mathematical statistics and in statistical 
physics, this is also known as detailed balance. 

the solutions to the corresponding deterministic equations (13.71) : yi{t) x and 
y2{t) X A^, respectively. 



7 Conclusion and discussion 

This paper has two intertwined threads. First, it reports a novel observation that for 
an SIS epidemic dynamics in a finite population with heterogeneous subgroups, 
there is an emergent inherent frequency, even though the corresponding determin- 
istic nonlinear dynamics for infinite population shows no indication whatsoever. 
This result provides new insight to the understanding of fluctuations in epidemio- 
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Figure 5: A typical stochastic realization of the SIS model. Here we considered two 
subgroups with k = 1,2. The parameters used in simulation are A^i = 50, = 50 (thus 
{k) = 1.5), and the transmission rate A = 1. The simulation is carried out according to 
Gillespie's algorithm ll48l . The thickened lines are orbits of the deterministic equations 
(13.71 ) with same A and Di = D2 = ^. 

logical data. 

Since the phenomenon is at the juncture between dynamics in terms of systems 
of deterministic ordinary differential equations (ODEs) and in terms of stochas- 
tic birth-death processes, the relation between these two types of mathematical 
models is rigorously investigated. We want to emphasize that in our approach, the 
nonlinear dynamics is an emergent, collective phenomenon of the stochastic pop- 
ulation dynamics. The second thread of the paper, therefore, is to advance a sys- 
tematic approach to nonlinear population dynamics based on individual's behavior 
with uncertainties, which can be represented in terms of probability distributions. 
From this starting point, the traditionally employed ODEs can be more explic- 
itly justified as the limiting behavior of infinitely large population for a nonlinear 
stochastic dynamics, called Delbriick-Gillespie process in cellular biochemistry 
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Il3l |46]|. Furthermore, a Gaussian like stochastic process can also be derived to 
account for the population fluctuations in the large, but finite, population. 

Indeed, on the level of individuals, epidemic processes are essentially stochas- 
tic, and the heterogeneity are inevitable. In the present study, we distinguish the 
stochastic, but statistically identical behavior (i.e., within each subgroup) from 
statistically non-identical behavior (i.e., between different subgroups). We show 
a more careful model building based "mechanisms" of the infection, though still 
quite crude, nevertheless can provide further insights into the dynamics on the 
population level. 

One might wonder why the circular dynamics disappears in deterministic sys- 
tem. The reason lies in the equality 9(^rn,n) = m+ln+i ■ have shown, the ODE 
dynamics equations are obtained from stochastic processes via the Law of Large 
Numbers, which is in order of A^^, where N is the total number of individuals in 
the population. In the limit of — )■ oo, 9(^m,n) 1- 

In summary, we proposed a microscopic, statistically heterogeneous contact 
process for individual-based epidemiological modeling. The fundamental stochas- 
tic events are the instantaneous "contacts" between two individuals and the recov- 
ery of an infectious individual. Within an infinitesimal time interval, these events 
are sequential and independent. The stochastic process is a Poisson flow. We 
applied this approach to two-state SIS dynamics. Following the method of Q ex- 
pansion, we obtained a system of deterministic ODEs for the densities for a large 
population, as well as a linear diffusion system near the stable steady state. 

There are two time scales in the dynamics of SIS epidemics. Compared with 
the time in which the infection reaches a "stationary pattern" among the sub- 
groups, the time for ultimate extinction is very long when A > Ac Il371 . We used 
the method of quasi-stationary approximation to study the long-time behavior on 
the "short time scale". 

Another assumption adopted here is the mass action law which requires well 
mixing. In population dynamic term, it is assumed that individuals are moving 
rapidly within a relatively small, highly "fluid" community. In a large spatial 
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scales, the contacts and spatial movements involve geographic factors; then a spa- 
tial model is required. In the applications of epidemic dynamic models to "virus 
infection" on the World Wide Web, however, this restriction is hardly there. 
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